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The Hamiltonian dynamics associated to classical, planar, Heisenberg XY models is investigated 
for two- and three-dimensional lattices. Besides the conventional signatures of phase transitions, 
here obtained through time averages of thermodynamical observables in place of ensemble averages, 
qualitatively new information is derived from the temperature dependence of Lyapunov exponents. 
A Riemannian geometrization of newtonian dynamics suggests to consider other observables of 
geometric meaning tightly related with the largest Lyapunov exponent. The numerical computation 
of these observables - unusual in the study of phase transitions - sheds a new light on the microscopic 
dynamical counterpart of thermodynamics also pointing to the existence of some major change 
in the geometry of the mechanical manifolds at the thermodynamical transition. Through the 
microcanonical definition of the entropy, a relationship between thermodynamics and the extrinsic 
geometry of the constant energy surfaces Eg of phase space can be naturally established. In this 
framework, an approximate formula is worked out, determining a highly non-trivial relationship 
between temperature and topology of the E_g. Whence it can be understood that the appearance 
of a phase transition must be tightly related to a suitable major topology change of the E_b- This 
contributes to the understanding of the origin of phase transitions in the microcanonical ensemble. 
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I. INTRODUCTION 



The present paper deals with the study of the microscopic Hamiltonian dynamical phenomenology associated to 
thermo dynamical phase transitions. This general subject is addressed in the special case of planar, classical Heisenberg 
XY models in two and three spatial dimensions. A preliminary presentation of some of the results and ideas contained 
in this paper has been already given in jjj. 

There are several reasons to tackle the Hamiltonian dynamical counterpart of phase transitions. On the one hand, 
we might wonder whether our knowledge of the already wide variety of dynamical properties of Hamiltonian systems 
can be furthcrly enriched by considering the dynamical signatures, if any, of phase transitions. On the other hand, 
it is a-priori conceivable that also the theoretical investigation of the phase transition phenomena could benefit of a 
direct investigation of the natural microscopic dynamics. In fact, from a very general point of view, we can argue 
that in those times where microscopic dynamics was completely unaccessible to any kind of investigation, statistical 
mechanics has been invented just to replace dynamics. During the last decades, the advent of powerful computers has 
made possible, to some extent, a direct access to microscopic dynamics through the so called molecular dynamical 
simulations of the statistical properties of "macroscopic" systems. 

Molecular dynamics can be either considered as a mere alternative to Monte Carlo methods in practical computa- 
tions, or it can be also seen as a possible link to concepts and methods (those of nonlinear Hamiltonian dynamics) 
that could deepen our insight about phase transitions. In fact, by construction, the ergodic invariant measure of the 
Monte Carlo stochastic dynamics, commonly used in numerical statistical mechanics, is the canonical Gibbs distribu- 
tion, whereas there is no general result that guarantees the ergodicity and mixing of natural (Hamiltonian) dynamics. 
Thus, the general interest for any contribution that helps in clarifying under what conditions equilibrium statistical 
mechanics correctly describes the average properties of a large collection of particles, safely replacing their microscopic 
dynamical description. 

Actually, as it has been already shown and confirmed by the results reported below, there are some intrinsically 
dynamical observables that clearly signal the existence of a phase transition. Notably, Lyapunov exponents appear 
as sensitive measurements for phase transitions. They are also probes of a hidden geometry of the dynamics, because 
Lyapunov exponents depend on the geometry of certain "mechanical manifolds" whose geodesic flows coincide with 
the natural motions. Therefore, a peculiar energy - or temperature - dependence of the largest Lyapunov exponent 
at a phase transition point also reflects some important change in the geometry of the mechanical manifolds. 

As we shall discuss throughout the present paper, also the topology of these manifolds has been discovered to play 
a relevant role in the phase transition phenomena (PTP). 

Another strong reason of interest for the Hamiltonian dynamical counterpart of PTP is related to the equivalence 
problem of statistical ensembles. Hamiltonian dynamics has its most natural and tight relationship with microcanon- 
ical ensemble. Now, the well known equivalence among all the statistical ensembles in the thermodynamic limit is 
valid in general in the absence of thermodynamic singularities, i.e. in the absence of phase transitions. This is not a 
difficulty for statistical mechanics as it might seem at first sight rather, this is a very interesting and intriguing 
point. 

The inequivalence of canonical and microcanonical ensembles in presence of a phase transition has been analytically 
shown for a particular model by Hertel and Thirring jj| , it is mainly revealed by the appearance of negative values of 
the specific heat and has been discussed by several authors |i|,|5| . 

The microcanonical description of phase transitions seems also to offer many advantages in tackling first order 
phase transitions ||, and seems considerably less affected by finite-size scaling effects with respect to the canonical 
ensemble description JtJ. This non-equivalence problem, together with certain advantages of the microcanonical 
ensemble, strenghtens the interest for the Hamiltonian dynamical counterpart of PTP. Let us briefly mention the 
existing contributions in the field. 

Butera and Caravati considering an XY model in two dimensions, found that the temperature dependence of the 
largest Lyapunov exponent changes just near the critical temperature T c of the Kosterlitz-Thouless phase transition. 
Other interesting aspects of the Hamiltonian dynamics of the XY model in two dimensions have been extensively 
considered in |9| , where a very rich phenomenology is reported. Recently, the behaviour of Lyapunov exponents has 
been studied in Hamiltonian dynamical systems: i) with long-range interactions flC|-|l^], ii) describing either clusters 
of particles or magnetic or gravitational models exhibiting phase transitions, in) in classical lattice field theories with 
O(l), 0(2) and 0(4) global symmetries in two and three space dimensions [OjM ], iv) in the XY model in two and 
three space dimensions [Q, v) in the "O - transition" of homopolymeric chains fig]. The pattern of A(T) close to the 
critical temperature T c is model-dependent. The behaviour o f Ly apunov exponents near the transition point has been 
considered also in the case of first- order phase transitions |lqjl7f| . It is also worth mentioning the very intriguing 
result of Ref. ]18||, where a glassy transition is accompanied by a sharp jump of A(T). 

X(T) always detects a phase transition and, even if its pattern close to the critical temperature T c is model- 
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dependent, it can be used as an order parameter - of dynamical origin - also in the absence of a standard order 
parameter (as in the case of the mentioned " 0-transition" of homopolymers and of the glassy transition in amorphous 
materials). This appears of great prospective interest also in the light of recently developed analytical methods to 
compute Lyapunov exponents (see Section IV). 

Among Hamiltonian models with long-range interactions exhibiting phase transitions, the most extensively studied 
is the mean-field XY model |pd| , ^9| - ^ll , whose equilibrium statistical mechanics is exactly described, in the thermo- 
dynamic limit, by mean-field theory]ll| . In this system, the theoretically predicted temperature dependence of the 
largest Lyapunov exponent A displays a non- analytic behavior at the phase transition point. 

The aims of the present paper are 

• to investigate the dynamical phenomenology of Kosterlitz-Thouless and of second order phase transitions in the 
2d and 3d classical Heisenberg XY models respectively; 

• to highlight the microscopic dynamical counterpart of phase transitions through the temperature dependence 
of the Lyapunov exponents, also providing some physical interpretation of abstract quantities involved in the 
geometric theory of chaos (in particular among vorticity, Lyapunov exponents and sectional curvatures of con- 
figuration space); 

• to discuss the hypothesis that phase transition phenomena could be originated by suitable changes in the topology 
of the constant energy hypersurfaces of phase space, therefore hinting to a mathematical characterization of phase 
transitions in the microcanonical ensemble. 

The paper is organized as follows: Sections II and 77/ are devoted to the dynamical investigation of the 2d and 3c? 
XY models respectively. In Section IV the geometric description of chaos is considered, with the analytic derivation 
of the temperature dependence of the largest Lyapunov exponent, the geometric signatures of a second-order phase 
transition and the topological hypothesis. Section V contains a presentation of the relationship between the extrinsic 
geometry and topology of the energy hypersurfaces of phase space and thermodynamics; the results of some numeric 
computations are also reported. Finally, Section VI is devoted to summarize the achievements reported in the present 
paper and to discuss their meaning. 



II. 2D XY MODEL 



We considered a system of planar, classical "spins" (in fact rotators) on a square lattice of N = n x n sites, and 
interacting through the ferromagnetic interaction V = —Ylu j\ JSi • Sj (where |S;| = 1). The addition of standard, 
i.e. quadratic, kinetic energy term leads to the following choice of the Hamiltonian 

H = | ~2~ + J t 2 ~ cos (ft+ij - <H,j) ~ cos (<li,j+i ~ lid)] | > W 

where j are the angles with respect to a fixed direction on the reference plane of the system. In the usual definition 
of the XY model both the kinetic term and the constant term 2JN are lacking; however, their contribution does 
not modify the thermodynamic averages (because they usually depend only on the configurational partition function, 
Zc = J rii=i dQi exp[— f3V(q)], the momenta being trivially integrable when the kinetic energy is quadratic). Thus, 
as we tackle classical systems, the choice of a quadratic kinetic energy term is natural because it corresponds to 
h y~lj—i |Si| 2 , written in terms of the momenta pij canonically conjugated to the lagrangian coordinates qij. The 
constant term 2JN is introduced to make the low energy expansion of Eq. (Q) coincident with the Hamiltonian of a 
system of weakly coupled harmonic oscillators. 

The theory predicts for this model a Kosterlitz-Thouless phase transition occurring at a critical temperature esti- 
mated around T c ~ J. Many Monte Carlo simulations of this model have been done in order to check the predictions 
of the theory. Among them, we quote those of Tobochnik and Chester [^] and of Gupta and Baillie [^3| which, on 
the basis of accurate numerical analysis, confirmed the predictions of the theory and fixed the critical temperature at 
T c = 0.89 (J = 1). 

The analysis of the present work is based on the numerical integration of the equations of motion derived from 
Hamiltonian ([!]) . The numerical integration is performed by means of a bilateral, third order, symplectic algorithm p4[ , 
and it is repeated at several values of the energy density e = E/N (E is the total energy of the system which depends 
upon the choice of the initial conditions). While the Monte Carlo simulations perform statistical averages in the 
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canonical ensemble, Hamiltonian dynamics has its statistical counterpart in the microcanonical ensemble. Statistical 
averages are here replaced by time averages of relevant observables. In this perspective, from the microcanonical 
definition of temperature l/T = dS/dE, where S is the entropy, two definitions of temperature are available: T = 
jj(K) (where K is the kinetic energy per degree of freedom), if S = log/H^i dqidpi<d(H(p, q) — E), where 0(-) 

is the Heaviside step function, and f = [(f - l) {K' 1 )]' 1 , if S = logJU.Zi d 9i d Pi S i H (P^) ~ E ) S T ( or T) 
are numerically determined by measuring the time average of the kinetic energy K per degree of freedom (or its 
inverse), i.e. T = lim^—nx, jj j f Q cLtK{t) (and similarly for T). There is no appreciable difference in the outcomes of 
the computations of temperature according to these two definitions. 

A. Dynamical analysis of thermodynamical observables 

1. Order parameter 

The order parameter for a system of planar "spins" whose Hamiltonian is invariant under the action of the group 
0(2), is the bidimensional vector 

n n n n 

M = ( E s "i> E S L-) = ( E E sin ^)' ( 2 ) 

which describes the mean spin orientation field. After the Mermin- Wagner theorem, we know that no symmetry- 
breaking transition can occur in one and two dimensional systems with a continuous symmetry and nearest-neighbour 
interactions. This means that, at any non-vanishing temperature, the statistical average of the total magnetization 
vector is necessarily zero in the thermodynamic limit. However, a vanishing magnetization is not necessarily expected 
when computed by means of Hamiltonian dynamics at finite N . In fact, statistical averages are equivalent to averages 
computed through suitable markovian Monte Carlo dynamics that a-priori can reach any region of phase space, 
whereas in principle a true ergodicity breaking is possible in the case of differentiable dynamics. Also an "effective" 
ergodicity breaking of differentiable dynamics is possible, when the relaxation times - of time to ensemble averages - 
are very fastly increasing with N |26| . 

This model has two integrable limits: coupled harmonic oscillators and free rotators, at low and high temperatures 
respectively. Hereafter, T is meant in units of the coupling constant J. 

For a lattice of N = 10 x 10 sites, Figure [j] shows that at low temperatures (T < 0.5)- being the system almost 
harmonic - we can observe a persistent memory of the total magnetization associated with the initial condition, which, 
on the typical time scales of our numeric simulations (10 6 units of proper time), looks almost frozen. 

By raising the temperature above a first threshold Tq ~ 0.6, the total magnetization vector - observed on the same 
time scale - starts rotating on the plane where it is confined. A further increase of the temperature induces a faster 
rotation of the magnetization vector together with a slight reduction of its average modulus. 

At temperatures slightly greater than 1, we observe that already at N — 10 x 10 a random variation of the direction 
and of the modulus of the vector M(i) sets in. 

At T > 1.2, we observe a fast relaxation and, at high temperatures (T ~ 10), a sort of saturation of chaos. 

At a first glance, the results reported in Fig. || could suggest the presence of a phase transition associated with 
the breaking of the 0(2) symmetry. In fact, having in mind the Landau theory, the ring-shaped distribution of the 
instantaneous magnetization shown by Fig. [j] is the typical signature of an 0(2)-broken symmetry phase and the 
spot-like patterns around zero are proper to the unbroken symmetry phase. 

The apparent contradiction of these results with the Mermin- Wagner theorem is resolved by checking whether the 
observed phenomenology is stable with N. Thus, some simulations have been performed at larger values of N. At 
any temperature, we found that the average modulus (|M(i)|) t of the vector M(t), computed along the trajectory, 
systematically decreases by increasing N. However, for temperatures lower than To, the iV-dependence of the order 
parameter is very weak, whereas, for temperatures greater than To, the iV-dependence of the order parameter is rather 
strong. In Fig. |2| two extreme cases (N = 10 x 10 and N = 200 x 200) are shown for T = 0.74. The systematic trend 
of (|M(t)|) toward smaller values at increasing TV is consistent with its expected vanishing in the limit N — > oo. 

At T = 1, Fig. shows that, when the lattice dimension is greater than 50 x 50, M(t) displays random variations 
both in direction (in the interval [0,27r]) and in modulus (between zero and a value which is smaller at larger N). 
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2. Specific heat 



By means of the recasting of a standard formula which relates the average fluctuations of a generic observable 
computed in canonical and microcanonical ensembles ]2^| , and by specializing it to the kinetic energy fluctuations, 
one obtains a microcanonical estimate of the canonical specific heat 



cv(T) = — 



k B d 

cv(e) = — 
T = T(e) 



1 - 



Nd (K 2 ) - (K) 2 ^' ' 



(3) 



where d is the number of degrees of freedom for each particle. Time averages of the kinetic energy fluctuations 
computed at any given value of the energy density e yield CV(T), according to its parametric definition in Eq.(||). 

^From the microcanonical definition 1/CV = dT(E) / dE of the constant volume specific heat, a formula can be 
worked out |25|, which is exact at any value of N (at variance with the expression ([§). It reads 



cv 



Cv 
N 



= {N-(N-2)(K)(K- 1 )Y 



(4) 



and it is the natural expression to be used in Hamiltonian dynamical simulations of finite systems. 

The numerical simulations of the Hamiltonian dynamics of the 2d XY model - computed with both Eqs.([|) and (^) 
- yield a cuspy pattern for cy(T) peaked at T ~ 1 (Fig. This is in good agreement with the outcomes of canonical 
Monte Carlo simulations reported in Ref. p2| , p3| , where a pronounced peak of cy(T) was detected at T ~ 1.02. 

By varying the lattice dimensions, the peak height remains constant, in agreement with the absence of a symmetry- 
breaking phase transition. 



3. Vorticity 



Another thermodynamic observable which can be studied is the vorticity of the system. Let us briefly recall that 
if the angular differences of nearby "spins" are small, we can suppose the existence of a continuum limit function 
0(r) that conveniently fits a given spatial configuration of the system. Spin waves correspond to regular patterns of 
0(r), whereas the appearance of a singularity in 6{r) corresponds to a topological defect, or a vortex, in the "spin" 
configuration. When such a defect is present, along any closed path C that contains the centre of the defect, one has 

<j)V6{r)-dr = 2nq, q = 0, ±1, ±2, . . . (5) 

indicating the presence of a vortex whose intensity is q. For a lattice model with periodic boundary conditions, there 
is an equal number of vortices and antivortices (i.e. vortices rotating in opposite directions). Thus, the vorticity of 
our model can be defined as the mean total number of equal sign vortices per unit volume. In order to compute 
the vorticity V as a function of temperature, we have averaged the number of positive vortices along the numerical 
phase space trajectories. On the lattice, r is replaced by the multi- index i and V M #.; = — eft, then the number 
of elementary vortices is counted: the discretized version of ^ V9 • dr = 1 amounts to one elementary vortex on a 
plaquette. Thus V is obtained by summing over all the plaquettes. 

Our results are in agreement with the values obtained by Tobochnik and Chester 
simulations with N = 60 x 60. 

As shown in Fig. [|, on the 10 x 10 lattice, the first vortex shows up at T ~ 0.6 and on the 40 x 40 lattice at 
T ~ 0.5, when the system changes its dynamical behavior, increasing its chaoticity (see next Subsection). At lower 
temperatures, vortices are less probable, due to the fact that the formation of vortex has a minimum energy cost. 
Below T ~ 1, the vortex density steeply grows with a power law V(T) ~ T 10 . The growth of V then slows down, until 
the saturation is reached at T ~ 10. 



22 1 by means of Monte Carlo 
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B. Lyapunov exponents and chaoticity 



The values of the largest Lyapunov exponent Ai have been computed using the standard tangent dynamics equations 
[see Eqs. ( |l0| ) and (|5^)], and are reported in Fig. [|. 

Below T ~ 0.6, the dynamical behavior is nearly the same as that of harmonic oscillators and the excitations of the 
system are only "spin- waves" . 

In the interval [0., 0.6], the observed temperature dependence Ai(T) ~ T 2 is equivalent to the Ai(e) ~ e 2 dependence 
(since at low temperature T(e) oc e), already found - analytic ally and numerically - in the quasi-harmonic regime of 
other systems and characteristic of weakly chaotic dynamics j28| . 

Above T ~ 0.6, vortices begin to form and correspondingly the largest Lyapunov exponent signals a "qualitative" 
change of the dynamics through a steeper increase vs. T. 

At T ~ 0.9, where the theory predicts a Kostcrlitz - Thoulcss phase transition, Ai(T) displays an inflection point. 

Finally, at high temperatures, the power law Ai(T) ~ T -1 / 6 is found. 



III. 3D XY MODEL 

In order to extend the dynamical investigation to the case of second-order phase transitions, we have studied a 
system described by an Hamiltonian having at the same time the main characteristics of the 2d model and the 
differences necessary to the appearance of a spontaneous symmetry-breaking below a certain critical temperature. 
The model we have chosen is such that the spin rotation is constrained on a plane and only the lattice dimension has 
been increased, in order to elude the "no go" conditions of the Mermin- Wagner theorem. This is simply achieved by 
tackling a system defined on a cubic lattice ofJV = nxnxn sites and described by the Hamiltonian 



H = y~] {^jp + J [3 - cos^i+i^fc - q itji k) 



cosfej+i,*; - q iij: k) - cosfai^fc+i - qi,j,k)]} ■ (6) 



A. Dynamical analysis of thermodynamical observables 

The basic thermodynamical phenomenology of a second-order phase transition is characterized by the existence of 
equilibrium configurations that make the order parameter bifurcating away from zero at some critical temperature T c 
and by a divergence of the specific heat cy(T) at the same T c . Therefore, this is the obvious starting point for the 
Hamiltonian dynamical approach. 



1. Order parameter 

Below a critical value of the temperature, the symmetry-breaking in a system invariant under the action of the 0(2) 
group, appears as the selection - by the average magnetization vector of Eq. (||)- of a preferred direction among all 
the possible, energetically equivalent choices. By increasing the lattice dimension, the symmetry breaking is therefore 
characterized by a sort of simultaneous "freezing" of the direction of the order parameter M and of the convergence 
of its modulus to a non-zero value. 

Figure ^ shows that in the 3d lattice, at T < 2, i.e. in the broken-symmetry phase (as we shall see in the following), 
the dynamical simulations yield a thinner spread of the longitudinal fluctuations by increasing N - that is, |M| 
oscillates by exhibiting a trend to converge to a non-zero value - and that the transverse fluctuations damp, "fixing" 
the direction of the oscillations. This direction depends on the initial conditions. 

Moreover, the dynamical analysis provides us with a better detail than a simple distinction between regular and 
chaotic dynamics. In fact, it is possible to distinguish between three different dynamical regimes (Fig. ||). 

At low temperatures, up to T ~ 0.8, one observes the persistency of the initial direction and of an equilibrium value 
of the modulus IMI close to one. 
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At 0.8 < T < 2.2, one observes transverse oscillations, whose amplitude increases with temperature. 

At T > 2.2, the order parameter exhibits the features typical of an unbroken symmetry phase. In fact, it displays 
fluctuations peaked at zero, whose dispersion decreases by increasing the temperature (bottom of Fig. ||) and, at a 
given temperature, by increasing the lattice volume (Fig. |9|a,b). 

We can give an estimate of the order parameter by evaluating the average of the modulus (|M(t)|) = p(T). At 
T < 2.2, the iV-dependence is given mainly by the rotation of the vector, while the longitudinal oscillations are 
moderate, as shown in Fig. [h]. At temperatures above T ~ 2.2, we observe the squeezing of p(T) to a small value. 

The existence of a second order phase transition can be recognized by comparing the temperature behavior and the 
A-dependence of the thermodynamic observables computed for the 2c? and the 3d models. Both systems exhibit the 
rotation of the magnetization vector and small fluctuations of its modulus when they are considered on small lattices. 
In the 2d model the average modulus of the order parameter is theoretically expected to vanish logarithmically with 
N, what seems qualitatively compatible with the weak N dependence shown in Fig. || whereas in the 3c? model we 
observe a stability with N of (|M|), suggesting the convergence to a non-zero value of the order parameter also in the 
limit N — ► oo, as shown in Fig. ^. 

T ~ 2.2 is an approximate value of the critical temperature T c of the second-order phase transition. This value will 
be refined in the following Subsection. No finite-size scaling analysis has been performed for two different reasons: i) 
our main concern is a qualitative phenomenological analysis of the Hamiltonian dynamics of phase transitions rather 
than a very accurate quantitative analysis, ii) finite-size effects are much weaker in the microcanonical ensemble than 
in the canonical ensemble M . 



2. Specific heat 

As in the 2d model, numerical simulations of the Hamiltonian dynamics have been performed with both Eqs.^) 
and (Eh. The outcomes show a cusplike pattern of the specific heat, whose peak makes possible a better determination 
of the critical temperature. By increasing the lattice dimension up to N — 15 x 15 x 15, the cusp becomes more 
pronounced, at variance with the case of the 2d model. Fig. O shows that this occurs at the temperature T c ~ 2.17. 



3. Vorticity 

The definition of the vorticity in the 3d case is not a simple extension of the 2c? case. Vortices are always defined 
on a plane and if all the "spins" could freely move in the three-dimensional space, the concept of vortices would 
be meaningless. For the 3c? planar (anisotropic) model considered here, vortices can be defined and studied on two- 
dimensional subspaces of the lattice. The variables Qij,k do not contain any information about the position of the 
plane where the reference direction to measure the angles qij,k is assigned. Dynamics is completely independent of 
this choice, which has no effect on the Hamiltonian. Moreover, as the Hamiltonian is symmetric with respect to the 
lattice axes, the three coordinate-planes are equivalent. This equivalence implies that vortices can contemporarily 
exist on three orthogonal planes. Though the usual pictorial representation of a vortex can hardly be maintained, its 
mathematical definition is the same as in the 2c? lattice case. Hence three vorticity functions exist and their average 
values - at a given temperature - should not differ, what is actually confirmed by numerical simulations. 

The vorticity function vs. temperature is plotted in Fig. [l2]. On a lattice of 10 x 10 x 10 spins, the first vortex is 
observed at T ~ 0.8. The growth of the average density of vortices is very fast up to the critical temperature, above 
which the saturation is reached. 



B. Lyapunov exponents and symmetry-breaking phase transition 

A quantitative analysis of the dynamical chaoticity is provided by the temperature dependence of the largest 
Lyapunov exponent. 

Figure |l3| shows the results of this computation. At low temperatures, in the limit of quasi-harmonic oscillators, 
the scaling law is again found to be X\(T) ~ T 2 and, at high temperatures, the scaling law is again \\{T) ~ y- 1 /^ 
as in the 2d case. In the temperature range intermediate between T ~ 0.8 and T c ~ 2.17, there is a linear growth 
of \i(T). At the critical temperature, the Lyapunov exponent exhibits an angular point. This makes a remarkable 
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difference between this system undergoing a second order phase transition and its 2d version, undergoing a Kosterlitz- 
Thouless transition. In fact, the analysis of the 2d model has shown a mild transition between the different regimes 



of Ai(T) (inset of Fig. 12), whereas in the 3d model this transition is sharper (inset of Fig. hj 



We have also computed the temperature dependence of the largest Lyapunov exponent of Markovian random 
processes which replace the true dynamics on the energy surfaces Y>e (see Appendix). The results are shown in Fig. 
[jjj . The dynamics is considered strongly chaotic in the temperature range where the patterns Ai(T) are the same for 
both random and differentiable dynamics, i.e. when differentiable dynamics mimics, to some extent, a random process. 
The dynamics is considered weakly chaotic when the value Ai resulting from random dynamics is larger than the value 
Ai resulting from differentiable dynamics. The transition from weak to strong chaos is quite abrupt. Figure |l4] shows 
that the pattern of the largest Lyapunov exponent computed by means of the random dynamics reproduces that of 
the true Lyapunov exponent at temperatures T > T c . This means that the setting in of strong thermodynamical 
disorder corresponds to the setting in of strong dynamical chaos. The "window" of strong chaoticity starts at T c and 
ends at T ~ 10. The existence of a second transition from strong to weak chaos is due to the existence, for T — > oo, 
of the second integrable limit (of free rotators), whence chaos cannot remain strong at any T > T c . 



IV. GEOMETRY OF DYNAMICS AND PHASE TRANSITIONS 



Let us briefly recall that the geometrization of the dynamics of ./V-degrees-of-freedom systems defined by a La- 
grangian C = K — V, in which the kinetic energy is quadratic in the velocities: K = ^aijtfcf , stems from the fact 
that the natural motions are the extrema of the Hamiltonian action functional Sh = J Cdt, or of the Maupertuis' 
action Sm — 2 J K dt. In fact, also the geodesies of Riemannian and pseudo-Riemannian manifolds are the extrema 
of a functional, the arc-length I = j ds, with ds 2 = gijdq l dq 3 . Hence, a suitable choice of the metric tensor allows 
for the identification of the arc-length with either Sh or Sm , and of the geodesies with the natural motions of the 
dynamical system. Starting from Sm, the "mechanical manifold" is the accessible configuration space endowed with 
the Jacobi metric p{J 

(g J ) ij = [E-V(q)}a ij , (7) 

where V(q) is the potential energy and E is the total energy. A description of the extrema of Hamilton's action Sh 
as geodesies of a "mechanical manifold" can be obtained using Eisenhart's metric |3(| on an enlarged configuration 
spacetime ({q° = t, q 1 , . . . , q N } plus one real coordinate q N+1 ), whose arc-length is 

ds 2 = -2V({q})(dq ) 2 + a i:j dq i dq j + 2dq°dq N+1 . (8) 

The manifold has a Lorentzian structure and the dynamical trajectories are those geodesies satisfying the condition 
ds 2 = Cdt 2 , where C is a positive constant. In the geometrical framework, the (in)stability of the trajectories is the 
(in)stability of the geodesies, and it is completely determined by the curvature properties of the underlying manifold 
according to the Jacobi equation [^9| , ^ll 



ds 2 3km ds ? ds 



R \ik m -T-Z -j- = , (9) 



whose solution £, usually called Jacobi or geodesic variation field, locally measures the distance between nearby 
geodesies; V/ds stands for the covariant derivative along a geodesic and R 1 - km are the components of the Riemann 
curvature tensor. Using the Eisenhart metric (||), the relevant part of the Jacobi equation (j^) is pS} ] 



2 c i 



d 2 £ 
dt 2 



R l ok oC = o , 2 = 1,...,^ (io) 



where the only non- vanishing components of the curvature tensor are Rmoj = d 2 V/ dqtdqj . Equation ( |l0| ) is the tangent 
dynamics equation, which is commonly used to measure Lyapunov exponents in standard Hamiltonian systems. Having 
recognized its geometric origin, it has been devised in Ref. p^j a geometric reasoning to derive from Eq. ( |l0| ) an effective 
scalar stability equation that, independently of the knowledge of dynamical trajectories, provides an average measure 
of their degree of instability. An intermediate step in this derivation yields 

^- + k R (t)?+6KW(t)? =0 , (11) 
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where Ur = Kr/N is the Ricci curvature along a geodesic defined as Kr = ^Rijcfqi, with v 2 = q l qi and R d = R k ik j, 

and 5K^ is the local deviation of sectional curvature from its average value . The sectional curvature is defined 

as = R ijkiCq j £ k q 1 / II £ l| 2 || 4 IP- 

Two simplifying assumptions are made: (i) the ambient manifold is almost isotropic, i.e. the components of the 
curvature tensor — that for an isotropic manifold (i.e. of constant curvature) are Rijkm = koidikgjm — gimgjk), 
k = const - can be approximated by Rijkm ~ k{t){gikgj m — gi m gjk) along a generic geodesic ~f(t); (ii) in the large 
TV limit, the "effective curvature" k(t) can be modeled by a gaussian and ^-correlated stochastic process. Hence, one 
derives an effective stability equation, independent of the dynamics and in the form of a stochastic oscillator equation 

S 

^ + [fc + wj(t)]V> = , (12) 

where ip 2 cx |£| 2 . The mean fc and variance Ok of k(t) are given by k — (Kr)/N and a 2 — ((Kr — (Kr)) 2 )/N, 
respectively, and the averages (•) are geometric averages, i.e. integrals computed on the mechanical manifold. These 
averages are directly related with microcanonical averages, as it will be seen at the end of Section V. r](t) is a gaussian 
(^-correlated random process of zero mean and unit variance. 

The main source of instability of the solutions of Eq. ( |l2|) , and therefore the main source of Hamiltonian chaos, is 
parametric resonance, which is activated by the variations of the Ricci curvature along the geodesies and which takes 
place also on positively curved manifolds [p2| . The dynamical instability can be enhanced if the geodesies encounter 
regions of negative sectional curvatures, such that kR + SK^ < 0, as it is evident from Eq. (pi]). 

In the case of Eisenhart metric, it is K R = AV = Etit 3 ^/ 5 ??) and r{2) = R OiOjC & / II £ || 2 = 
(d 2 V/dq t dq 3 )£, z £, j /||£|| 2 . The exponential growth rate A of the quantity ip 2 + ip 2 of the solutions of Eq. (|l2"l), is 
therefore an estimate of the largest Lyapunov exponent that can be analytically computed. The final result reads |28] 



A 2fc . l 2 /64fc3 



A= \2a 2 T + \ -Jl+Aafr 2 , (13) 



where r = ny/k^ / (2^ ko(k + a k ) + Trer*.); in the limit ak/k <C 1 one finds A cx a 2 . 

A. Signatures of phase transitions from geometrization of dynamics 

In the geometric picture, chaos is mainly originated by the parametric instability activated by the fluctuating 
curvature felt by geodesies, i.e. the fluctuations of the (effective) curvature are the source of the instability of 
the dynamics. On the other hand, as it is witnessed by the derivation of Eq. ([l2]) and by the equation itself, a 
statistical-mechanical-like treatment of the average degree of chaoticity is made possible by the geometrization of 
the dynamics. The relevant curvature properties of the mechanical manifolds are computed, at the formal level, 
as statistical averages, like other thermodynamic observables. Thus, we can expect that some precise relationship 
may exist between geometric, dynamic and thermodynamic quantities. Moreover, this implies that phase transitions 
should correspond to peculiar effects in the geometric observables. 

In the particular case of the 2d XY model, the microcanonical average kinetic energy (K) and the average Ricci 
curvature (Kr) computed with the Eisenhart metric are linked by the equation 

I N d 2 v\ N 

Kr =\Y1 )= 2J Y, ( C0S fe+iJ - *j) + cos(q id+1 - q ld )) = 2( J -(V)) , (14) 



1,3 = 1 J 



so that 



H = Ne=(K) + (V)^-^ = e-2J+^- 1 ^- . (15) 

Being the temperature defined as T = 2(K)/N (with ks = 1) and being d = 1 (because each spin has only one 
rotational degree of freedom), from Eq.(0) it follows that 



If lal/Ny 1 
2\ 2 T 2 J 



cv = - 2 [i-- n ^r) ■ (16) 



9 



In the special case of these XY systems, it is possible to link the specific heat and the Ricci curvature by inserting 
Eq.(|l5|) into the usual expression for the specific heat at constant volume. Thus, one obtains the equation 

1 d(K R )(T) 

CV = ~2N dT ■ (17) 

The appearance of a peak in the specific heat function at the critical temperature has to correspond to a suitable 
temperature dependence of the Ricci curvature. 

In the 3d model, the potential energy and the Ricci curvature are proportional, according to: jt{V) =3 — -^(Kr) . 

Another interesting point is the relation between a geometric observable and the vorticity function in both models. 
As already seen in previous sections, the vorticity function is a useful signature of the dynamical chaoticity of the 
system. From the geometrical point of view, the enhancement of the instability of the dynamics with respect to 
the parametric instability due to curvature fluctuations, is linked to the probability of obtaining negative sectional 
curvatures along the geodesies (as discussed for Id XY model in Ref. (2^]). In fact, when vortices are present in the 
system, there will surely be two neighbouring spins with an orientation difference greater than 7r/2, such that, if i,j 
and i + are their coordinates on the lattice, it follows that 

Qi+1,3 - %j > 2 cos W+lJ ~ %j) < ' ( 18 ) 

The sectional curvature relative to the plane defined by the velocity v along a geodesic and a generic vector £ _L v is 

k<»= y m ^ d9) 



i,j,k,l=X 



For the 2d XY model, it is 
J 



J N 

K(2) = T^W E { C0S (%+U - qiM i+1 ' j - + oo S (q iJ+1 - qi,j)[C' j+1 - C' 3 ?} ■ (20) 



Thus, a large probability of having a negative value of the cosine of the difference among the directions of two close 
spins corresponds to a larger probability of obtaining negative values of the sectional curvatures along the geodesies; 
here for £ the geodesic separation vector of Eq. ([h]) is chosen. 

In the 3d model, the sectional curvature relative to the plane defined by the velocity v and a generic vector £ _L v 

is 

J N 



Kir- , 



+ cosfe J+1 , fe - g^kHf - e j 'T + cos(<h Afc+x - qi, jlk )[e' J ' k+ - e ,3 'T} ( 2l ) 

and again the probability of finding negative values of K^ 2 ' along a trajectory is limited to the probability of finding 
vortices. 

The mean values of the geometric quantities entering Eq. (|l2|) can be numerically computed by means of Monte 
Carlo simulations or by means of time averages along the dynamical trajectories. In fact, due to the lack of an explicit 
expression for the canonical partition function of the system, these averages are not analytically computable. For 
sufficiently high temperatures, the potential energy becomes negligible with respect to the kinetic energy, and each 
spin is free to move independently from the others. Thus, in the limit of high temperatures, one can estimate the 
configurational partition function Zq — J*" Hi dq\e~^ v ^ by means of the expression 

« JV N 

Zc = e~ 2tiJN I dq itj exp{/3J ^ [cosfe+ij - q itj ) + cos(q lJ+1 - qij)]} 

i,3 = l i,j=l 
,7T n N 

^ e -2f)jN I j-j duijdvi i3 ■ exp{/3J ^2 [ cos ( u i,j) + cos(vi,j)]} (22) 

i,j=l i,j=l 

after the introduction of Uij = q^ij — qij and Vij = qi t j+i—qij as independent variables. In this way, some analytical 
estimates of the average Ricci curvature ko(T) and of its r.m.s. fluctuations <J 2 {T) have been obtained for the 2d 
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model (Fig. |l5|). For temperatures above the temperature of the Kosterlitz-Thouless transition, these estimates are 
in agreement with the numerical computations onaiV= lOx 10 lattice. It is confirmed that Hamiltonian dynamical 
simulations, already on rather small lattices, are useful to predict, with a good approximation, the thermodynamic 
limit behavior of relevant observables. Moreover, the good quality of the high temperature estimate gives a further 
information: at the transition temperature, the correlations among the different degrees of freedom are destroyed, 
confirming the strong chaoticity of the dynamics. 

The same high temperature estimates of ko(T) and cr^(T) have been performed for the 3d system. In Fig. [l|, 
the numerical determination of a\ (T) shows the appearance of a very pronounced peak at the phase transition point 
which is not predicted by the analytic estimate, whereas the average Ricci curvature ko(T) is in agreement with the 
analytic values of the high temperature estimate, computed by spin decoupling, above the critical temperature, as in 
the 2d model. 



B. Geometric observables and Lyapunov exponents 



We have seen that the largest Lyapunov exponent is sensitive to the phase transition and at the same time we 
know that it is also related to the average curvature properties of the "mechanical manifolds" . Thus, the geometric 
observables ko(T) and cr^(T) above considered can be used to estimate the Lyapunov exponents, as well as to detect 
the phase transition. 



In principle, by means of Eq. ( 13 ) , one can evaluate the largest Lyapunov exponent without any need of dynamics, but 
simply using global geometric quantities of the manifold associated to the physical system. For 2d and 3d XY models, 
fully analytic computations are possible only in the limiting cases of high and low temperatures. Microcanonical 
averages of ko and <j\ at arbitrary T have been numerically computed through time averages. We can call this hybrid 
method semi- analytic. 

In Fig. [l7[ the results of the semi-analytic prediction of the Lyapunov exponents for the 2d model are plotted vs. 
temperature and compared with the numerical outcomes of the tangent dynamics. As one can see, the prediction 



formulated on the basis of Eq.(13) underestimates the numerical values given by the tangent dynamics. The semi- 
analytic prediction can be improved by observing that the replacement of the sectional curvature fluctuation SK^ in 
Eq.(|Tl|) with a fraction of the Ricci curvature [which underlies the derivation of Eq.(|l2])] underestimates the frequency 
of occurrence of negative sectional curvatures, which was already the case of the Id XY model p8|| . The correction 
procedure can be implemented by evaluating the probability P(T) of obtaining a negative value of the sectional 
curvature along a generic trajectory and then by operating the substitution 

1 + P(l )a 

The parameter a is a free parameter to be empirically estimated. Its value ranges from 100 to 200, without appreciable 
differences in the final result. It resumes the non trivial information about the more pronounced tendency of the 
trajectories towards negative sectional curvatures with respect to the predictions of the geometric model describing 
the chaoticity of the dynamics. 

The probability P(T) is estimated through the occurrence along a trajectory of negative values of the sum of the 
coefficients that appear in the definition of [Eqs.@ and @)] 

, , /-^ Q(~ cos(g fc+ i,; - g k ,i) - cos(g fc ,; + i - q K i)) exp[-/?y(q)] []|^ dq k ,i 

exp[-/3V(q)] n fe ,i=i dq k ,l 

averaged over all the sites Vk,l <E (1, . . . , N); 9 is the step function. 

Alternatively, owing to the already remarked relation between vorticity and sectional curvature K^ 2 \ P(T) can be 
replaced by the average density of vortices 

« r >-iwr < 25 > 

where a a free parameter. Actually, in the 2d model, the two corrections, one given by Eq.([2^) with P(T) of Eq. 
(pi|), the other given by Eq.(^5|) with the vorticity function in place of P(T), convey the same information. The 
semi-analytic predictions of Ai(T) with correction are reported in Fig. [l^. 
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In the limits of high and low temperatures, Ai(T) can be given the analytic forms Ai(T) ~ J 1 — 1/6 a t high tempera- 
ture, and Ai(T) ~ T 2 at low temperature. In the former case, the high temperature approximation (|2|) is used, and in 
the latter case the quasi-harmonic oscillators approximation is done. The deviation of Ai(T) from the quasi-harmonic 
scaling, starting at T ~ 0.6 and already observed to correspond to the appearance of vortices, finds here a simple 
explanation through the geometry of dynamics: vortices are associated with negative sectional curvatures, enhancing 
chaos. 

By increasing the spatial dimension of the system, it becomes more and more difficult to accurately estimate the 
probability of obtaining negative sectional curvatures. The assumption that the occurrence of negative values of the 
cosine of the difference between the directions of two nearby spins is nearly equal to P(T), is less effective in the 3c? 
model than in the 2d one. Again, the vorticity function can be assumed as an estimate of P(T) [Eq. (p5|)]. The 
quality of the results has a weak dependence upon the parameter a. The correction remains good, with a belonging 
to a broad interval of values (100 -i- 200). In the limits of high and low temperatures, the model predicts correctly the 
same scaling laws of the 2d system. 

In Fig. 03 the semi-analytic predictions for the Lyapunov exponents, with and without correction, are plotted vs. 
temperature together with the numerical results of the tangent dynamics. It is noticeable that the prediction of Eq. 
© is able to give the correct asymptotic behavior of the Lyapunov exponents also at low temperatures, the most 
difficult part to obtain by means of dynamical simulations. 



C. A topological hypothesis 



We have seen in Fig. that a sharp peak of the Ricci-curvature fluctuations <y 2 (T) is found for the 3d model 
in correspondence of the second order phase transition, whereas, for the 2d model, o- 2 (T) appears regular and in 
agreement with the theoretically predicted smooth pattern. On the basis of heuristic arguments, in Refs. 0JT^] we 
suggested that the peak of a 2 observed for the 3d XY model, as well as for 2d and 3d scalar and vector lattice ip 4 
models, might originate in some change of the topology of the mechanical manifolds. In fact, in abstract mathematical 
models, consisting of families of surfaces undergoing a topology change - i.e. a loss of diffeomorphicity among them 
at some critical value of a parameter labelling the members of the family, we have actually observed the appearance 
of cusps of a 2 K at the transition point between two subfamilies of surfaces of different topology, K being the Gauss 
curvature. 

Actually, for the mean- field XY model, where both o- 2 (T) and Ai(T) have theoretically been shown to loose analytic- 
ity at the phase transition point, a direct evidence of a "special" change of the topology of equipotential hypersurfaces 
of configuration space has been given (3^] . Other indirect and direct evidences of the actual involvement of topology 
in the deep origin of phase transitions have been recently given |}4 

In the following Section we consider the extension of this topological point of view about phase transitions from 
equipotential hypersurfaces of configuration space to constant energy hypersurfaces of phase space. 



35 1 for the lattice tp 4 model. 



V. PHASE SPACE GEOMETRY AND THERMODYNAMICS. 



In the preceding Section we have used some elements of intrinsic differential geometry of submanifolds of configu- 
ration space to describe the average degree of dynamical instability (measured by the largest Lyapunov exponent). 
In the present Section we are interested in the relationship between the extrinsic geometry of the constant energy 
hypersurfaces E^ and thermodynamics. 

Hereafter, phase space is considered as an even-dimensional subset T of M. 2N and the hypersurfaces Eg = 
{(pi, . . . ,pn-, Qi, ■ ■ ■ , Qn) £ M\H {pi, . . . ,PNi Qi, ■ ■ ■ , In) = E} are manifolds that can be equipped with the standard 
Riemannian metric induced from M. 2N . If, for example, a surface is parametrically defined through the equations x i — 
x l (z 1 , . . . , z k ), i — 1, . . . , 2N, then the metric gij induced on the surface is given by gijiz 1 , . . . , z k ) — Yln=i lfTTfl~- 
The geodesic flow associated with the metric induced on Eg from M. 2N has nothing to do with the Hamiltonian flow 
that belongs to S^. Nevertheless, it exists an intrinsic Riemannian metric gs of phase space T such that the geodesic 
flow of gs, restricted to Eg, coincides with the Hamiltonian flow (gs is the so called Sasaki lift to the tangent bundle 
of configuration space of the Jacobi metric gj that we mentioned in a preceding Section) . 

The link between extrinsic geometry of the Eg and thermodynamics is estabilishcd through the microcanonical 
definition of entropy 
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t 2N 



is the metric induced from 



where da — y/ det(g)dxi...dx2N-i is the invariant volume element of T,e C 
and X1...X2N-1 are the coordinates on Te- 

Let us briefly recall some necessary definitions and concepts that are needed in the study of hypersurfaces of 
euclidean spaces. 

A standard way to investigate the geometry of an hypersurface S™ 1 is to study the way in which it curves around 
in R m+1 : this is measured by the way the normal direction changes as we move from point to point on the surface. 
The rate of change of the normal direction N at a point x € £ is described by the shape operator L x (y) = — V V N = 
— (ViVi • v, . . . , VAT m+ i ■ v), where v is a tangent vector at x and V v is the directional derivative of the unit normal N. 
As L x is an operator of the tangent space at x into itself, there are m independent eigenvalues p6| «i(x), . . . , n m (x), 
which are called the principal curvatures of X at x. Their product is the Gauss- Kronecker curvature: Kq(x) = 
YliLi K i(x) = det(i x ), and their sum is the so-called mean curvature: M\{x) — — Y^Li K%{x). The quadratic form 
L x (y) ■ v, associated with the shape operator at a point x, is called the second fundamental form of X at x. 

It can be shown [pl[ that the mean curvature of the energy hypersurfaces is given by 



Mi (a;) 



1 



2N- 1 



VH(x) 
\VH(x)\ 



(27) 



where VH(x)/||VH(x)|| is the unit normal to Eg 
(d/dpi, . . . , d/dq n), whence the explicit expression 



at a given point x = (pi, . . . ,Pn, ?!)•••; Qn), an( i V 



(2A^~ 1) Mi = - 



+ 



1 



\VH\ 
1 

Wm 



d 2 v 



/ d 2 V 



\dqidqjj \dqi 



dV 



dV 

dm 



(28) 



where i, j are multi-indices according to the number of spatial dimensions. 

The link between geometry and physics stems from the microcanonical definition of the temperature 



1 



dS 



i <m v {E) 



T dE 9, V (E) dE ' 
where we used Eq.@ with k B = 1, v = 2N - 1, and ft u (E) = J da/\\VH\\. /From the formula (37 

d k 



dE k 



j a da) (E') = j A k (a)da 



(29) 



(30) 



where a is an integrable function and A is the operator A(a) 
result 



1 
T 



1 dVL v {E) _ 1 



Q V {E) dE 



da 



NH 



\VH\ 



AH 



\\YH\ 



, it is possible to work out the 



VH\ 



1 



da 



Ml 



NH\\ \\VH\\ 



(31) 



where M\ 



V(Vi//|| VH || ) is directly proportional to the mean curvature (p7|). In the last term of Eq.(^TJ) we 
have neglected a contribution which vanishes as 0(X/N). Eq. ( pl| ) provides the fundamental link between extrinsic 
geometry and thermodynamics |p8| . In fact, the microcanonical average of M*/||VH||, which is a quantity tightly 
related with the mean curvature of He, gives the inverse of the temperature, whence other important thermodynamic 
observables can be derived. For example, the constant volume specific heat 



1 _ dT{E) 
CV ~ dE 



using Eq.(E9j), yields 



C v 



dS_ 
~d~E 



fd 2 S\ 
\dE 2 ) 



(32) 



(33) 
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becoming at large N 

C v = 





r 


1 


d f 


da / 


\l|Vff|| 


I mc 






l|VF|| ( 



Wh\\ 



R(E) 



)-(M 



(34) 



where the subscript mc stands for microcanonical average, and R(E) stands for the quantities of order 0(1/ N) 
neglected in the last term of Eq.(p9|) (a-priori, its derivative can be non negligible and has to be taken into account). 
Eq. (|34| ) highlights a more elaborated link between geometry and thermodynamics: the specific heat depends upon 
the microcanonical average of M* /||V_ff|| and upon the energy variation rate of the surface integral of this quantity. 

Remarkably the relationship between curvature properties of the constant energy surfaces £ e and thermodynamic 
observables given by Eqs.(^9|) and (|34j) can be extended to embrace also a deeper and very interesting relationship 
between thermodynamics and topology of the constant energy surfaces. Such a relationship can be discovered through 
a reasoning which, though approximate, is highly non-trivial, for it makes use of a deep theorem due to Chern and 
Lashof |p9[. As ||Vff|| = {J2iPi + ['ViT^(q')] 2 } 1 / 2 is a positive quantity increasing with the energy, we can write 



1 

T 



1 dtt v 
~^~d~E 



1 



da Mt 



D < E) <k 



da M ± , 



(35) 



where we have introduced the factor function D(E) in order to extract the total mean curvature J" s da Mi; D(E) has 
been numerically found to be smooth and very close to (1/|| ViJ|| 2 ) mc (see Section V A and Fig. |l9|). Then, recalling 
the expression of a multinomial expansion 



E 



mi 



! . . 



and identifying the Xi with the principal curvatures ki, one obtains 

V 

Ml = v\ Y[ h + R = v\K + 

where K = Y[iki is the Gauss-Kronecker curvature, and R is the sum 
(rife = 1, Vfc). Using v\ ~ v v e~ v \f 7 ha>, 



Ml ~ v v eT v ^^NK + R 
is obtained. The above mentioned theorem of Chern and Lashof states that 



■■^■■■K" , (36) 

R , (37) 
without the term with the largest coefficient 

(38) 



f \K\d0>V<A[®Z\^bi{V E 

i=Q 



(39) 



i.e. the total absolute Gauss-Kronecker curvature of a hypersurface is related with the sum of all its Betti numbers 



6i(Ss). The Betti numbers are diffeomorphism invariants of fundamental topological meaning 
sum is also a topologic invariant. §^ is a hypersphere of unit radius. Combining Eqs. 
on £e, we obtain 



101, therefore their 



\Ml\da 



j| ) and (39j) and integrating 
\K\da+ I \R\da>A^b l {Y JE )+TZ{E) , (40) 



with the shorthands A = v v e~ v V ol(Sl) and K = J^ E \R\da. 

Now, with the aid of the inequality J Wfl^^dfi > \\ J fd/iW 1 ^'" , we can write 



f \Mt\da = { \MZ\V v dtr> f Ml da 



l/v 



(41) 



If Mi > everywhere on then 

largely prevailing pi]], J Sb Ml da 
and therefore 



I: 

l/v 



,, Ml da 



l/v 



SE <Ml\da 

l/v 



l/v 



whence, in the hypothesis that Mi > is 



J s |Af 1 / |dcr] . Under the same assumption, J E Mida ~ J s |Mi|dcr 



14 



Mi da 



Finally, 



1 M V \ X I V dn > 




l/u , 




1 M v dn 




[ 1 M v 1 dn\ 
J \ m l 1 aa J 










1 dfl v 


1 M\ \ 




f da M* 


~ Q. v dE 


\IIVtf||/ m< 




k B \\VH\\ \\VH\\ 



> 



Aj2Uz E ) + n{E) 



i=0 



l/u 



(42) 



D(E) 



9. 



da Mi 



v Js 



> 



D(E) 



Aj2 b *(z E ) + n(E) 



l/u 



(43) 



Equation ( |43| ) has the remarkable property of relating the microcanonical definition of temperature of Eq. ( |35| ) with a 
topologic invariant of S^;. The Betti numbers can be exponentially large with TV [for example, in the case of TV-tori 



T^, they are bk — \T)]> so that the quantity (J] b^) 1 ^ can converge, at arbitrarily large TV, to a non-trivial limit (i.e. 
different from one). Thus, even though the energy dependence of 1Z is unknown, the energy variation of &»(5je) must 
be mirrored - at any arbitrary TV - by the energy variation of the temperature. By considering Eq. (p4|) in the light 
of Eq.(f43|), we can expect that some suitably abrupt and major change in the topology of the can reflect into the 
appearance of a peak of the specific heat, as a consequence of the associated energy dependence of bk(^E) and of its 
derivative with respect to E. In other words, we see that a link must exist between thermodynamical phase transitions 
and suitable topology changes of the constant energy submanifolds of the phase space of microscopic variables. The 
arguments given above, though in a still rough formulation, provide a first attempt to make a connection between 
the topological aspects of the microcanonical description of phase transitions and the already proposed topological 
hypothesis about topology changes in configuration space and phase transitions |l],[l4|,|33]-[55| . 

Direct support to the topological hypothesis has been given by the analytic study of a mean-field XY model |33| ] and 
by the numerical computation of the Euler characteristic x °f the equipotential hypersurfaces X„ of the configuration 
space in a 2d lattice cp 4 model |35| . The Euler characteristic is the alternate sum of all the Betti numbers of a manifold, 
so it is another topological invariant, but it identically vanishes for odd dimensional manifolds, like the Se. In Ref. p5[ , 
x(£„) neatly reveals the symmetry-breaking phase transition through a sudden change of its variation rate with the 
potential energy density v. A sudden "second order" variation of the topology of the H v appears in both Refs. [|33|j35[ 
as the requisite for the appearance of a phase transition. These results strenghten the arguments given in the present 
Section about the role of the topology of the constant energy hypersurfaces. In fact, the larger is N, the smaller are the 
relative fluctuations (S 2 V) 1 ^ 2 / (V) and (S 2 !^) 1 / 2 /{K) of the potential and kinetic energies respectively. At very large 



N, the product manifold X 



N-l 



x § 



N-l 



with v = (V) and t = (K), v + 1 = E, is a good model manifold to represent 



the part of that is overwhelmingly sampled by the dynamics and that therefore constitutes the effective support 



The kinetic energy submanifolds E> t = {(pi, . . . ,Pn) € 



= *} 



of the microcanonical measure on E^. 
are hyperspheres. 

In other words, at very large TV the microcanonical measure mathematically extends over a whole energy surface but, 
as far as physics is concerned, a non-negligible contribution to the microcanonical measure is in practice given only by 
a small subset of an energy surface. This subset can be reasonably modeled by the product manifold x S^ -1 , 

because the total kinetic and total potential energies - having arbitrarily small fluctuations, provided that TV is large 
enough - can be considered almost constant. Thus, since § t _1 at any t is always an hypersphere, a change in the 
topology of directly entails a change of the topology of x Sf^ 1 , that is of the effective model-manifold 

for the subset of ~Ee where the dynamics mainly "lives" at a given energy E. 

At small TV, the model with a single product manifold is no longer good and should be replaced by the non-countable 



union 



-ijV-l 



eici 



xS 



N-l 
E—v- 



with v assuming all the possible values in a real interval T. From this fact the smoothing 



of the energy dependence of thermodynamic variables follows. Nevertheless, the geometric and topologic signals of the 
phase transition can remain much sharper than the thermodynamic signals also at small TV (< 100), as it is witnessed 
by the 2d lattice tp 4 model [H||]. 

Finally, let us comment about the relationship between intrinsic geometry, in terms of which we discussed the 
geometrization of the dynamics, and extrinsic geometry, dealt with in the present Section. 

The most direct and intriguing link is estabilished by the expression for microcanonical averages of generic observ- 
ablcs of the kind A(q), with q = (qi, . . . , qjy), 



{A ) mc — 



1 



^2N(E) JH(p,q)<E 



d N pd N qA{q) = 



1 



Vol{M E ) J V (q)<E 



d N q[E-V{q)\ N l 2 A{q) , 



(44) 



where M E = {q e R N \V(q) < E}. Eq. @ is obtained by means of a Laplace-transform method |25| ; it is remarkable 
that [E — V(q)] N = det(gj), where gj is the Jacobi metric whose geodesic flow coincides with newtonian dynamics (see 
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Section V), therefore d N q[E — V(q)] N ^ 2 = d N q \J det(gj) is the invariant Riemannian volume element of [Me,9j)- 
Thus, 

v ,L s I d N q [E V(q)] N / 2 A(q) = ] [ d N q y/deJ(g^) A(q) , (45) 

Vol(M E ) J v(q) < E Vol{M E ) J Me 

which means that the microcanonical averages (A(q) } mc can be expressed as Riemannian integrals on the mechanical 
manifold (Me,<?j). 

In particular, this also applies to the microcanonical definition of entropy 

S = k B log / d N pd N q = k B log ( E dE' [ , (46) 



H(p,q)<E 



which is alternative to that given in Eq.(26), though equivalent to it in the large N limit. We have 



S = k B log 



= fe s log/ d N q ^Jdet(gj) + const . , (47) 

JMe 

where the last term gives the entropy as the logarithm of the Riemannian volume of the manifold. 

The topology changes of the surfaces that are to be associated with phase transitions, will deeply affect 

also the geometry of the mechanical manifolds (M B ,gj) and (M x R 2 ,^) and, consequently, they will affect the 
average instability properties of their geodesic flows. In fact, Eq.(13) links some curvature averages of these manifolds 
with the numeric value of the largest Lyapunov exponent. Loosely speakin g, m ajor topology changes of will 
affect microcanonical averages of geometric quantities computed through Eq. ( |44| ) , likewise entropy, computed through 
Eq.@. 

Thus, the peculiar temperature patterns displayed by the largest Lyapunov exponent at a second-order phase 
transition point - in the present paper reported for the 3c? XY model, in Ref. [Q reported for lattice <p 4 models 
- appear as reasonable consequences of the deep variations of the topology of the equipotential hypersurfaces of 
configuration space. 

We notice that topology seems to provide a common ground to the roots of microscopic dynamics and of thermo- 
dynamics and, notably, it can account for major qualitative changes simultaneously occurring in both dynamics and 
thermodynamics when a phase transition is present. 



A. Some preliminary numerical computations 



Let us briefly report on some preliminary numerical computations concerning the extrinsic geometry of the hyper- 
surfaces He in the case of the 3d XY model. 

The first point about extrinsic geometry that we numerically addressed was to check whether the inverse of the 
temperature, that appears in Eq.(p5|), can be reasonably factorized into the product of a smooth "deformation factor" 
D(E) and of the total mean curvature J E M\da. To this purpose, the two independently computed quantities 
(l/||ViJ|| 2 ) mc and D{E) = [J^Jda/WV H\\)(M* /\\V H\\)] / [J^ daA'h] are compared in Fig. [Oj, showing that actually 
J SE (da/\\VH\\)(Ml/\\VH\\) ~ (1/|| V£T|| 2 ) mc / Se doM x . In other words, D(E) ~ (1/|| ViJ|| 2 ) mc and no "singular" 
feature in its energy pattern seems to exist, what suggests that J s da Mi has to convey all the information relevant 
to the detection of the phase transition. There is no reason to think that the validity of the factorization given in 
Eq.(|35|) is limited to the special case of the XY model. 

The other point that we tackled concerns an indirect quantification of how a phase space trajectory curves around 
and knots on the to which it belongs. We can expect that the way in which an hypersurface is "filled" by 
a phase space trajectory living on it will be affected by the geometry and the topology of the ~Ee- In particular, we 
computed the normalized autocorrelation function of the time series Mi[x(f)] of the mean curvature at the points of 
YIe visited by the phase space trajectory, that is, the quantity 

T(T) = (SM 1 (t + T)SM 1 (t)) t , (48) 
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where 8M\(t) — Mi(t) — (Mi(t')) t > is the fluctuation with respect to the average (the "process" M\{t) is supposed 
stationary). Our aim was to highlight the extrinsic geometric-dynamical counterpart of a symmetry-breaking phase 
transition. 

The practical computation of T(r) proceeds by working out the Fourier power spectrum |Mi(w)| 2 of Afi[a;(t)], 
obtained by averaging 15 spectra computed by an FFT algorithm with a mesh of 2 15 points and a sampling time 
At = 0.1. Some typical results for T(t), obtained at different temperatures, are reported in Fig.^o| The patterns T(t) 
display a first regime of very fast decay, which is not surprising because of the chaoticity of the trajectories at any 
energy, followed by a longer tail of slower decay. An autocorrelation time r corr can be defined through the first intercept 
of T(r) with an almost-zero level (r = 0.01). In Fig.^l] we report the values of r corr so defined vs. temperature. In 
correspondence of the phase transition (whose critical temperature is marked by a vertical dotted line), r corr changes 
its temperature dependence: by lowering the temperature, below the transition T corr (T) rapidly increases, whereas 
it mildly decreases above the transition. Below T ~ 0.9, where the vortices disappear, the autocorrelation functions 
of Mi look quite different and it seems no longer possible to coherently define a correlation time. This result has 
an intuitive meaning and confirms that the phase transition corresponds to a change in the microscopic dynamics, 
as already signaled by the largest Lyapunov exponent; however, notice that the correlation times T corr (T) are much 
longer than the inverse values of the corresponding Ai(T). Qualitatively, Ai(T) and T~^. r (T) look similar, however 
the two functions are not simply related. 



VI. DISCUSSION AND PERSPECTIVES 

The microscopic Hamiltonian dynamics of the classical Heisenberg XY model in two and three spatial dimensions has 
been numerically investigated. This has been possible after the addition to the Heisenberg potentials of a standard 
(quadratic) kinetic energy term. Special emphasis has been given to the study of the dynamical counterpart of 
phase transitions, detected through the time averages of conventional thermodynamic observables, and to the new 
mathematical concepts that are brought about by Hamiltonian dynamics. 

The motivations of the present study are given in the Introduction. Let us now summarize what are the outcomes 
of our investigations and comment about their meaning. There are three main topics, tightly related one to the other: 

• the phenomenological description of phase transitions through the natural, microscopic dynamics in place of the 
usual Monte Carlo stochastic dynamics; 

• the investigation, in presence of phase transitions, of certain aspects of the (intrinsic) geometry of the mechanical 
manifolds where the natural dynamics is represented as a geodesic flow; 

• the discussion of the relationship between the (extrinsic) geometry of constant energy hypersurfaces of phase 
space and thermodynamics. 

About the first point, we have found that microscopic Hamiltonian dynamics very clearly evidences the presence of 
a second order phase transition through the time averages of conventional thermodynamic observables. Moreover, 
the familiar sharpening effects, at increasing N, of the specific heat peak and of the order parameter bifurcation are 
observed. The evolution of the order parameter with respect to the physical time (instead of the fictitious Monte 
Carlo time) is also accessible, showing the appearance of Goldstone modes and that, in presence of a second order 
phase transition, there is a clear tendency to the freezing of transverse fluctuations of the order parameter when N is 
increased. The "freezing" is observed together with a reduction of the longitudinal fluctuations, i.e. the rotation of 
the magnetization vector slows down, preparing the breaking of the 0(2) symmetry at N — ► oo. At variance, when a 
Kosterlitz-Thouless transition is present, at increasing TV the magnetization vector has a faster rotation and a smaller 
norm, preparing the absence of symmetry-breaking in the N — > oo limit as expected. 

Remarkably, to detect phase transitions, microscopic Hamiltonian dynamics provides us with additional observ- 
ables of purely dynamical nature, i.e. without statistical counterpart: Lyapunov exponents. Similarly to what we and 
other authors already reported for other models (see Introduction), also in the case of the 3d XY model a peculiar 
temperature pattern of the largest Lyapunov exponent shows up in presence of the second order phase transition, 
signaled by a "cuspy" point. By comparing the patterns Ai(T) given by Hamiltonian dynamics and by a suitably 
defined random dynamics respectively, we suggest that the transition between thermodynamically ordered and dis- 
ordered phases has its microscopic dynamical counterpart in a transition between weak and strong chaos. Though 
a — posteriori physically reasonable, this result is far from obvious, because the largest Lyapunov exponent measures 
the average local instability of the dynamics, which a — priori has little to do with a collective, and therefore global, 
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phenomenon such as a phase transition. The effort to understand the reason of such a sensitivity of Ai to a second 
order phase transition and to other kinds of transitions, as mentioned in the Introduction, is far reaching. 

Here we arrive to the second point listed above. In the framework of a Riemannian geometrization of Hamiltonian 
dynamics, the largest Lyapunov exponent is related to the curvature properties of suitable submanifolds of configura- 
tion space whose geodesies coincide with the natural motions. In the mathematical light of this geometrization of the 
dynamics, and after the numerical evidence of a sharp peak of curvature fluctuations at the phase transition point, 
the peculiar pattern of Ai(T) is due to some major change occurring to the geometry of mechanical manifolds at the 
phase transition. Elsewhere, we have conjectured that indeed some major change in the topology of configuration 
space submanifolds should be the very source of the mentioned major change of geometry. 

Thus, we have made a first attempt to provide an analytic argument supporting this topological hypothesis (third 
point of the above list). This is based on the appearance of a non trivial relationship between the geometry of constant 
energy hypersurfaces of phase space with their topology and with the microcanonical definition of thermodynamics. 
Even still in a preliminary formulation, our reasoning already seems to indicate the topology of energy hypersurfaces 
as the best candidate to explain the deep origin of the dynamical signature of phase transitions detected through 
Ai(T). ^ 

The circumstance, mentioned in the preceding Section, of the persistence at small N of geometric and topologic 
signals of the phase transition that are much sharper than the thermodynamic signals is of prospective interest for 
the study of phase transition phenomena in finite, small systems, a topic of growing interest thanks to the modern 
developments - mainly experimental - about the physics of nuclear, atomic and molecular clusters, of conformational 
phase transitions in homopolymers and proteins, of mesoscopic systems, of soft-matter systems of biological interest. 
In fact, some unambiguous information for small systems - even about the existence itself of a phase transition - could 
be better obtained by means of concepts and mathematical tools outlined here and in the quoted papers. Here we 
also join the very interesting line of thought of Gross and collaborators about the microcanonical description 

of phase transitions in finite systems. 

Let us conclude with a speculative comment about another possible direction of investigation related with this 
signature of phase transitions through Ly apunov exponents. In a field-theoretic framework, based on a path-integral 
formulation of classical mechanics [p^-[f5[ , Lyapunov exponents are defined through the expectation values of suitable 
operators. In the field-theoretic framework, ergodicity breaking appears related to a supersymmetry breaking [fj3f , 
and Lyapunov exponents are related to mathematical objects that have many analogies with topological concepts 
& 

The new mathematical concepts and methods, that the Hamiltonian dynamical approach brings about, could hope- 
fully be useful also in the study of more "exotic" transition phenomena than those tackled in the present work. Besides 
the above mentioned soft-matter systems, this could be the case of transition phenomena occurring in amorphous and 
disordered materials. 
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VIII. APPENDIX A 

Let us briefly explain how a random markovian dynamics is constructed on a given constant energy hypersurface 
of phase space. The goal is to compare the energy dependence of the largest Lyapunov exponent computed for 
the Hamiltonian flow and for a suitable random walk respectively. One has to devise an algorithm to generate a 
random walk on a given energy hypersurface such that, once the time interval At separating two successive steps is 
assigned, the average increments of the coordinates are equal to the average increments of the same coordinates for 
the differentiable dynamics integrated with a time step At. In other words, the random walk has to roughly mimick 
the differentiable dynamics with the exception of its possible time-correlations. 

One starts with a random initial configuration of the coordinates qi, i = 1, 2, . . . , N, uniformly distributed in the 
interval [0,27r], and with a random gaussian-distributed choice of the coordinates pt. The random pseudo-trajectory 
is generated according to the simple scheme 
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(Pi)(fe+i)At ^ (Pi)fcAt + a p G t , k At , 

where At is the time interval associated to one step k <— > k + 1 in the markovian chain, G^j 
random numbers with zero expectation value and unit variance; the parameters a q and otj 



(49) 

are gaussian distributed 
are the variances of the 



processes (qi)k and (pi)fc. These variances are functions of the energy per degree of freedom e. They have to be set 
equal to the numerically computed average increments of the coordinates obtained along the differentiable trajectories 
integrated with the same time step At, that is 




(?i(*H 


- At)- 
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" 1 w " 
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n 1/2 



(50) 



where T is the temperature. Then, in order to make minimum the energy fluctuations around any given value of the 
total energy, a criterium to accept or reject a new step along the markovian chain has to be assigned. A similar problem 
has been considered by Creutz, who developed a Monte Carlo microcanonical algorithm p6fl , where a "Maxwellian 
demon" gives a part of its energy to the system to let it move to a new configuration, or gains energy from the system, 
if the new proposed configuration produces an energy lowering. If the demon does not have enough energy to allow 
an energy increasing update of the coordinates, no coordinate change is performed. In this way, the total energy 
remains almost constant with only small fluctuations. As usual in Monte Carlo simulations, it is appropriate to fix 
the parameters so as the acceptance rate of the proposed updates of the configurations is in the range 30% - 60%. 

A reliability check of the so defined random walk, and of the adequacy of the phase space sampling through the 
number of steps adopted in each run, is obtained by computing the averages of typical thermodynamic observables of 
known temperature dependences. 

An improvement to the above described "demon" algorithm has been obtained through a simple reprojection on 
Tie of the updated configurations |^9| ; the coordinates generated by means of (^) are corrected with the formulae 



qi(kAt) h-> qi{kAt) + 



(ff)A£ 



(51) 



x R (kAt) 



p l {kAt) i-> p l (kAt) 



PiAE 



x R (kAt) 



where AE is the difference between the energy of the new configuration and the reference energy, and xu(kAt) denotes 
the random phase space trajectory. At each assigned energy, the computation of the largest Lyapunov exponent Af 
of this random trajectory is obtained by means of the standard definition 



1 " 

X? = lim — — log 



||C((fc + l)A*)|| 



.. a A! ^ " " ||C(fcAt)|| 



(52) 



where ((t) = £(t)) is given by the discretized version of the tangent dynamics 



&((fc + l)At) - 2&(fcAf) + U(k - l)At) 
At 2 



8 2 V 

dqidq-j 



&(*At)=0 



(53) 



x a (kAt) 



For wide variations of the parameters (At and acceptance rate), the resulting values of are in very good 
agreement. Moreover, the algorithm is sufficiently stable and the final value of Af is independent of the choice of the 
initial condition. 

A more refined algorithm could be implemented by constructing a random markovian process qitu) = 
[qi(tk): ■ ■ • i <7iv(ifc)] performing an importance sampling of the measure d/i — [E — V (q)] N dq in configuration space. 
In fact, similarly to what is reported in Eq.Q), one has [||| f H , p d N pd N q — const J V ^ <E d N q [E— V(q)\ N / 2 ~ l . 
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A random process obtained by sampling such a measure - with the additional property of a relation between the av- 
erage increment and the physical time step At as discussed above - would enter into Eq.(|5^) to yield Xf. However, 
this would result in much heavier numerical computations (with some additional technical difficulty at large N) which 
was not worth in view of the principal aims of the present work. 
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M.Cerruti-Sola, C.Clementi, M.Pettini 
Hamiltonian dynamics.... Fig. 1 
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FIG. 1. The magnetization vector M(t) computed along a trajectory for the 2d XY model at different temperatures on a 
lattice of N = 10 x 10. Each point represents a vector M(t) at a certain time t. 
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Hamiltonian dynamics Fig- 2 
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FIG. 2. The magnetization vector M(t) at the temperature T = 0.74, corresponding to the specific energy e = 0.8 and 
computed in a time interval At = 10 J , with a random initial configuration, on lattices of N = 10 x 10 (external points) and of 
N = 200 x 200 (internal points). 
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FIG. 3. The magnetization vector M(t) at the temperature T = 1, corresponding to the energy e = 1.2, computed in a time 
interval At = 10 5 , with a random initial configuration on lattices of a) N = 10 x 10, b) N = 50 x 50, c) N = 100 x 100 and 
iV = 200 x 200 sites, respectively. 
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FIG. 4. Specific heat at constant volume computed by means of Eq. (^) on a lattice of N = 10 x 10 (open circles) and 
of TV = 15 x 15 (full triangles). Starlike squares refer to specific heat values computed by means of Eq. (J3|) on a lattice of 
N = 10 x 10 . 
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FIG. 5. Vorticity function (plotted in a) linear scale and b) logarithmic scale) computed at different temperatures for 
lattices of N = 10 X 10 (open circles) and N — 40 x 40 (full circles). The results of the Monte Carlo simulations for a lattice 
of N — 60 x 60 (crosses) are from Kfl. The dashed line represents the power law V(t) ~ T 10 . 
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FIG. 6. The largest Lyapunov exponents computed on different lattice sizes: N = 10 x 10 (starred squares), N — 20 x 20 
(open triangles), N = 40 x 40 (open stars), iV = 50 x 50 (open squares) and N = 100 x 100 (open circles). In the inset, symbols 
have the same meaning. 
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Hamiltonian dynamics.... Fig- 7 
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FIG. 7. The magnetization vector M(£), computed at the temperature T = 1.7, on lattices of different sizes. By increasing 
the lattice dimensions, the longitudinal fluctuations decrease. The time interval At = 3.5 x 10 4 — 8 x 10 4 is the same for the 
four simulations. 
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FIG. 8. The magnetization vector M(t) computed at different temperatures on a lattice of iV = 10 x 10 x 10 spins. 
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FIG. 9. The magnetization vector M(t) computed at the temperature T = 2.22 (slightly higher than the critical value) on 
lattices of a) N = 10 x 10 x 10 and b)N = 15 x 15 x 15, respectively. The time interval At = 0.5 x 10 4 — 1.5 x 10 4 is the same 
for both simulations. 
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FIG. 10. The dynamical order parameter, defined as the average of the modulus |M(t)| along a trajectory, computed on 
lattices of N = 10 x 10 x 10 (full circles) and N = 15 x 15 x 15 (open circles). 
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FIG. 11. Specific heat at constant volume for the 3d model, computed by means of Eq. (^) on lattices of7V = 8x8x8 
(open triangles), N — 10 x 10 x 10 (open circles), N = 12 x 12 x 12 (open stars) and N = 15 x 15 x 15 (open squares) . Full 
circles refer to specific heat values computed by means of Eq. (|^) on a lattice of N = 10 x 10 x 10. The dashed line points out 
the critical temperature T c ~ 2.17 at which the phase transition occurs. 
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FIG. 12. Vorticity function at different temperatures along a dynamical trajectory on a lattice of N = 10 x 10 x 10 sites. 
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FIG. 13. The largest Lyapunov exponents computed at different temperatures for the 3d model. Numerical results are for 
lattices of iV = 10 x 10 x 10 (open circles) and iV = 15 x 15 x 15 (open stars). In the inset, symbols have the same meaning. 
The dashed line points out the temperature T c ~ 2.17 of the phase transition. The solid line puts in evidence the departure of 
Ai (T) from the linear growth. 
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FIG. 14. The largest Lyapunov exponents computed by means of the random dynamics algorithm (full circles) are 
i comparison with those computed by means of the standard dynamics (open stars) for a lattice of N — 10 x 10 x 10. 
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FIG. 15. Time average of Ricci curvature (open circles) and its r.m.s. fluctuations (full circles) at different temperatures 
computed for a lattice of N = 40 x 40 sites. Solid lines are the analytic estimates obtained from a high temperature expansion. 
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FIG. 16. Time average of Ricci curvature (open triangles) and its r.m.s. fluctuations (full triangles) computed at different 
temperatures for a lattice of N — 10 x 10 x 10. Open circles and full diamonds refer to a lattice size of N — 15 x 15 x 15. Solid 
lines are the analytic estimates in the limit of high temperatures. The dashed line points out the temperature T c ~ 2.17 of the 
phase transition. 
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FIG. 17. Analytic Lyapunov exponents computed for the 2d model by means of Eq.(|13[) without correction (dots) and 
incorporating the correction that accounts for the probability of obtaining negative sectional curvatures (full squares) for a 
lattice size of N = 40 x 40 are plotted in comparison with the numerical values of Fig. ^| The dashed lines are the asymptotic 
behaviors at high and low temperatures in the thermodynamic limit. 
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FIG. 18. Analytic Lyapunov exponents computed for the 3d model by means of Eq.(|13[) without correction (dots) and 
incorporating the correction that accounts for the probability of obtaining negative sectional curvatures (full circles) are plotted 
in comparison with the numerical values of Fig. [13j. The dashed lines are the asymptotic behaviors at high and low temperatures 
in the thermodynamic limit. 
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FIG. 19. The deformation factor D{E) = [/ (da/\\ VH\\)(Mf/\\ VH\\)] / [/ E doM{\ of Eq. (||) (open circles) is plotted 
vs. energy density E/N and compared to the quantity (1/||VJ/|| 2 ) (open triangles). N = 10 x 10 x 10. 
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FIG. 20. The normalized autocorrelation functions T(r) are plotted vs. time r for a lattice of N = 10 x 10 x 10 and for 
four different values of the temperature (from top to bottom: T = 0.49, 1.28, 1.75, 2.16). 
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FIG. 21. Autocorrelation times r corr are plotted vs. temperature T. The vertical dashed line points out the temperature 
T c ~ 2.17 at which the phase transition occurs. 
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